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ABSTRACT 

We have developed a technique to map the three-dimensional structure of the 
local interstellar medium using a maximum entropy reconstruction technique. A 
set of column densities N to stars of known distance can in principle be used to 
recover a three-dimensional density field n, since the two quantities are related 
by simple geometry through the equation N = C • n, where C is a matrix charac- 
terizing the stellar spatial distribution. In practice, however, there is an infinte 
number of solutions to this equation. We use a maximum entropy reconstruc- 
tion algorithm to find the density field containing the least information which 
is consistent with the observations. The solution obtained with this technique 
is, in some sense, the model containing the minimum structure. We apply the 
algorithm to several simulated data sets to demonstrate its feasibility and success 
at recovering "real" density contrasts. 

This technique can be applied to any set of column densities whose end points 
are specified. In a subsequent paper we shall describe the application of this 
method to a set of stellar color excesses to derive a map of the dust distribution, 
and to soft X-ray absorption columns to hot stars to derive a map of the total 
density of the interstellar medium. 

Subject headings: ISM: general 



1. Introduction 



Given a set of projected density measurements of columns with known end points, one 
can, at least in principle, recover the three dimensional distribution of absorbing material. 
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One way to accomplish this deprojection is to simply draw isocolumn contours (Sfeir et al. 
1999; Welsh et al. 1999). Here we propose an alternative technique - maximum entropy 
reconstruction (mer) - and apply a computer implementation of the algorithm to simulated 
data sets. We test its ability to recover structure in the interstellar medium (ism) by vary- 
ing several parameters and the functional form of the entropy. We then briefly sketch its 
application to a number of archived data sets, to be presented in a subsequent paper. 



The column or projected density of absorbing material iV(x) along the line of sight 
between an observer and a star at position r* is defined as 



where n(x) is the number density of absorbers at position x, usually measured in cm -3 . It 
is n(x), the map of the ISM, that we seek. Although the inversion of equation 1 is formally 
trivial, we know only N/r+, and so there is an infinite number of functions n(x) which 
satisfy the observed iV(x*). We therefore take a discrete approach, approximating the ISM 
as a coarse-binned density field. 

Consider a set of data consisting of iV* stars of known distance, to which accurate 
intervening column densities have been measured. We write the set of column densities to 
these N* stars as the vector N = (Aq, N 2 , . . . N^). Let us divide the space occupied by 
these stars into N c cells of size s. (In practice the cells need not be cubic - more generally 
we have s x , s y and s z . For example, it is sometimes useful to set s z to a large number to 
build a planar model of the ISM.) The model we seek consists of the set of number densities 
of these cells (n, a vector of length N c ) which will reproduce the observed column densities. 
The observations form a constraint equation 



where C is a matrix with N c rows and columns describing the line-of-sight/cell intersec- 
tions for each star. 

For the moment, let us assume that the observer is located at the corner shared by 8 
cells. (In the actual implementation of this algorithm we place the observer at the center of 
a cell, however, to enhance the numerical stability of several matrix inversions.) If s is small 



2. Maximum Entropy Technique 




(1) 



N = C-n, 



(2) 
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enough, no two lines of sight will intersect the same cell, and there will be an infinite number 
of solutions to this equation. It is clear, however, that some solutions will be preferred over 
others. For example, let us imagine that our data set consists of only one stellar position and 
intervening column density. It is easy to imagine two different models - the minimum and 
maximum structure configurations (Figure 1) - which reproduce the observations. Figure la 
shows the material spread out evenly along the line of sight, whereas in Figure lb the material 
inhabits only one cell. In the absence of any additional information both models reproduce 
the observations equally well. Information theory tells us that the former is to be preferred 
to the latter, however, since it makes the fewest assumptions about the underlying structure 
- i.e. it contains less information. Since our goal is to uncover structure in the ISM, the 
most conservative approach is to find the solution containing the minimum structure. Thus 
we seek a configuration which is a solution to equation 2 and which contains the minimum 
amount of information. 

The quantity of information contained in a map is usually measured using the Shannon 
entropy (Shannon 1948; Kapur & Kesavan 1992). To construct the map with the fewest 
assumptions about the underlying structure, but consistent with a set of observations, we 
invoke Jaynes' principle of maximum entropy (Jaynes 1957). The nomenclature "entropy" 
was adopted because the Shannon entropy was shown by Jaynes (1957) to be identical to 
the classical Boltzmann entropy. 

In general, a map can be characterized by a distribution function / = /(t, x, v), where 
/ is the mass per unit volume per velocity interval. The amount of information contained 
in the distribution is characterized by the Boltzmann H function, 



H = J /log/rfxrfv (3) 
which is related to the thermodynamic entropy S through 



Here V is the volume covered by the map, and k is Boltzmann's constant. 

The choice of grid resolution is fundamental to the character of the reconstructed map. 
We will choose a cell size s such that i>ism • ^obs <C s, where vism is the typical velocity of 
material in the ISM and t Q b s is the timescale spanned by the observations. Thus we can 
ignore any time dependence in the phase space density, making it a function of cell position 
only. In addition, vism <C c throughout the Galactic ISM, and so the absorptive properties 
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of the material in each cell are independent of the velocity of the material. Thus we replace 
/(x, v, t) with the coarse-binned pj. Scaling the entropy by the mean particle mass and 
Boltzmann's constant we obtain 



where rii is the number density of absorbers in cell i. The analogy between rti and probability 
Pi of finding a particle in cell i leads immediately to the identification of equation 5 as the 
Shannon entropy. We shall, however, refer to this form as the "Boltzmann entropy" of the 
map, in keeping with past astrophysical applications of the maximum entropy principle. 

In the present problem we are not restricted to using the Boltzmann entropy. Boltz- 
mann's H theorem was originally formulated to show that the thermodynamic entropy is 
the only non-decreasing function of the time-dependent distribution function for an ideal 
gas (Boltzmann 1896; Tolman 1938). In collisionless systems, however, there is an infi- 
nite number of H functions which exhibits this behavior (Tremaine, Henon, & Lynden-Bell 
1986). These are of the form 



where C is any convex function of the distribution function. Following Richstone & Tremaine 
(1988) we shall use the term "entropy" to refer to any function of this form, although, in 
the context of information, "uncertainty" might be a more appropriate term, as Kapur & 
Kesavan (1992) point out. 

The ISM is an extremely complex physical system. Although it can be considered a 
"collisional" system in some sense, astrophysical processes such as radiative cooling, energy 
injection by supernovae, and self-regulating star formation render this characterization of 
its self-interaction too simplistic. Thus there are no compelling arguments based upon first 
principles which favor one form of the entropy over another. We shall instead turn to 
numerical experimentation and convenience to guide our choice of the entropy. We shall 
return to this topic in Section 3.1. 

Regardless of the form of entropy we ultimately adopt, we seek seek the solution of 
equation 2 where S is a maximum, or, equivalently, we seek the map containing the minimum 
information. (Although there are many configurations wherein S is minimized, there is only 
one, for a given grid, where it is maximized.) This problem is readily solved using the method 




(5) 



i=i 




(6) 



i=i 
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of Lagrange multipliers (Jaynes 1957). For example, Richstone & Tremaine (1988) adapted 
this technique to solve the problem in a stellar dynamical context. We adopt this approach 
here, recasting it for the problem at hand. 

The (generally nonlinear) variational equation we must solve is 



SS — 5N • L = (7) 

where the vector L = (L 1; L 2 , . . . , L Nc ) is the set of Lagrange multipliers to be determined 
later. In the limit of small variations we have 



VS - (VN)* • L = (8) 

where the operator V is defined as 

^ ~ 9 „ d d 

V = ex— + e 2 — + ... + e N - 9 

orii on 2 ohn c 

and (VN)* is the transpose of VN. Since the column density is related to the number density 
through equation 2, VN = C, and we have then 

V,S-C*-L = (10) 

To solve this equation we guess a solution no and iterate using the Newton-Raphson 
method. Expanding about the initial guess and keeping terms only to first order yields 

VSU = VS\ no + V(VS| no ) • 5n (11) 
Inserting equation 11 into equation 10, and solving for 5n yields 

5n = (OSU- 1 • C* • L — (^L,)- 1 • VS\ no (12) 

where <$S = V(VS'). To make use of this equation we must first calculate the Lagrange 
multipliers. Dotting both sides with C and solving for L we obtain 

L = (C ■ OSI- 1 • Cr 1 • SN + (C • ^U, 1 ■ C')- 1 • C ■ ^I^ 1 • VS\ no (13) 
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Here we have identified C • Sn as 8N, the difference between the observations and the 
model column densities corresponding to the guessed solution n . Equation 12 is applied 
iteratively, with the Lagrange multipliers calculated anew at each step using equation 13, 
until a satisfactory solution is obtained. The computational expense is dominated by the 
inversion of (}S, an N c x N c matrix. 

The density of material in each cell is physically (though not mathematically) restricted 
to non-negative values. In order to impose this condition upon our solution we require that 
each component of Sn, as calculated in the code by equation 12, be no greater in magnitude 
than the current value of density in that cell. This guarantees that each cell has a non- 
negative density at every step in the iterative reconstruction and in the final maximum 
entropy solution. 

3. Simulated Data Sets 

When we apply the mer algorithm to a real data set, we will not know whether the 
resultant map is close to the actual ISM structure. (We can be sure that the reconstructed 
map is very different from the true structure if the iteration scheme fails to converge, how- 
ever.) To ascertain how well the algorithm reconstructs underlying distribution of absorbing 
matter we must test its performance upon simulated data sets. We postulate a 3D density 
structure n(x) for the ISM and a stellar spatial distribution, and calculate N(x+) for each 
star using equation 1. We then apply the mer algorithm to the faked data set N(x*) and 
compare the maximum entropy configuration with the simulated density structure. 

There are three aspects of the scheme which need to be explored using these simulated 
data. These are (1) the form of the entropy used in the maximization routine, (2) the size 
of the spatial cells used, for a given sampling density, and (3) the sampling density of the 
stars in the data set, compared with the scale length of the smallest structures in the ISM. 

3.1. Functional Form of the Entropy 

A strict interpretation of information theory dictates that the entropy we should max- 
imize is the Shannon entropy (i.e. the classical Boltzmann entropy), in the form given by 
equation 5. However, in the spirit of experimentation, we also ran tests for three other forms 
of the entropy as well. Two of these are "true" entropies in the the sense of Tremaine, Henon, 
& Lynden-Bell (1986), following the form of equation 6. They are 
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(14) 



i=i 



and 



(15) 



i=l 



The former, because it is a simple quadratic, should converge in only one iteration in the 
Newton-Raphson scheme. The fourth function we tried was 



This function is not of the form given by equation 6, and so we refer to it as a "pseudo- 
entropy". Although it is not of the form discussed by Tolman (1938) and Tremaine, Henon, 
& Lynden-Bell (1986), like these other forms it possesses a single global maximum - a 
property which, as it turns out, is probably sufficient to guarantee its utility. The four 
entropy types are compared in Figure 2. 

The simulated density field used for this experiment is shown in Figure 3. It consists of 
an ambient medium with a density of ~ 0.1 cm -3 , with two structures superimposed upon 
it. One is a wall of material with a gaussian amplitude ~ 150 cm~ 3 ; the other is a small 
cloud with gaussian amplitude ~ 70 cm~ 3 . The entire map has dimensions of 1017.5 pc x 
632.5 pc. 

Throughout this density field we distribute 458 stars in a random fashion (but with a 
minimum separation) and use the calculated columns to try to reconstruct the field. MER 
maps for the four entropy types given by equations 5, 14, 15, and 16 are shown in Figure 4. 
Each reconstruction uses a cell size of 27.5 pc x 27.5 pc, for a total grid size of 37 x 23. Only 
806 of the 851 cells are active in the reconstruction, however, because many edge cells are 
not pierced by a stellar line of sight. (These cells, shown in dark gray, are obvious toward 
the edges because they contain no stars.) We have set a column match tolerance of 10 for 
convergence termination, although in practice this will be set by the column measurement 
accuracy. (It should be noted that the cell size was optimized for this example.) Table 1 
lists the mean column error, rms column error, the mean density error, and the rms density 
error, for each of the reconstructions. Quantities with the mer subscript obviously refer to 
properties of the reconstruction. Column densities with no subscript represent the input data, 
while number densities with no subscript represent the density of the "actual" configuration 
at the cell center. 




(16) 



i=l 
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Two conclusions can be drawn from this experiment. The first is that mer is an effective 
way to recover real density contrasts without introducing significant spurious structure. Both 
structures are clearly reproduced by the algorithm with the appropriate amplitudes. And in 
no case do we find extraneous structure on scales larger than a pixel. 

The second observation is that the mer algorithm is quite insensitive to the form of the 
entropy used. The reconstructions created using the four different entropy types are nearly 
identical; in fact, one would be hard-pressed to distinguish one from another without referring 
to Table 1. Even the pixel-sized variations occur in the same location, suggesting that the 
"noise" in each reconstructed map is due to the peculiarities of the stellar distribution with 
respect to the cell grid, rather than tied to a specific form of the entropy. We thefore 
conclude that it is best to use the quadratic entropy, to save on computational expense. If 
we normalize the time to calculate the entropy and its derivatives to the that of the quadratic 
form, we have t Bo itz = 1-3, t cxp o = 1-5, and t pscu do = 2.0. 

In practice, program run-times are actually much shorter using the quadratic entropy 
because the mer algorithm usually converges in one or two iterations, compared with four 
to seven for the other entropy types. These run-time savings are greatest when the grid 
contains a large number of cells (and by extention a large number of stars, since they scale 
linearly in optimized configurations). In the 458-star simulation, a single iteration of the 
code running on a 450 MHz Pentium II processor took 52 s when using a grid of 684 cells, 
and 1230 s with a 2285-cell grid (see Section 3.2). So while using the quadratic entropy may 
only save a few minutes for a small number of stars, the difference can be several hours for 
a large data set. 

3.2. Cell size 

We next sought to understand the importance of cell size for a given stellar sampling 
density. The number density of the ISM should really be calculated differentially from the 
column data through equation 1; the MER algorithm, however, uses a differencing scheme. 
Since the number density of absorbers in a particular cell is constant across the cell in 
the reconstructed map, we expect there to be numerical instabilities when the cell size is 
comparable to or larger than the size scale of structures in the underlying ISM. If the cell is 
so large that it contains two stars, and a non- negligible density gradient, there may actually 
be no solution to the constraint equation, regardless of any entropy considerations. 

We examined the effect of cell size modulation on the reconstructed map, for a fairly 
dense simulated stellar distribution. For this set of experiments we use the simulated density 
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map shown in Figure 5. It contains three clouds of various size scale and maximum density. 
From largest to smallest, the gaussian amplitudes are ~ 100, 40, and 35 cm -3 . The entire 
map has dimensions of 1012.5 pc x 607.5 pc, with 666 stars distributed throughout. We 
varied the cell size from 27.5 pc on a side, down to 12.5 pc, in increments of —2.5 pc. These 
mer maps are shown in Figure 6. 

Table 2 gives the error statistics for each reconstruction. We omit the 27.5 pc cell mer 
map because the reconstruction failed - the algorithm was terminated after one iteration 
because the solution was divergent. The reason for this is that several cells contain two 
stars and a sufficiently steep density gradient, resulting in zero solutions of the constraint 
equation. Each of the reconstructions shown in Figure 6 fit the data reasonably well (the 
column tolerance was set to 10~ 3 ). As the cell size decreases, however, the solution begins 
to show radial striations as more and more cells are intersected by only one or two stellar 
lines of sight. In the limiting case of infinitesimal cells, every active cell in the grid is pierced 
only once; in essence, each individual line of sight reduces to the configuration illustrated 
in Figure 1. Thus material is smeared evenly along each absorption path, and the N* x N c 
optimization problem reduces to N+ separate ctkN c problems (52k=i ak = -0- Since no cell is 
intersected by more than one line of sight, these paths no longer "communicate" with each 
other, and their individual entropy contributions decouple during the optimization. 

For this particular case the solution is optimized (in the sense of matching the "true" 3D 
density) for a cell/star ratio of about 7/4. However, barring peculiarities in stellar position 
within each cell, (N mer — N) — > as the cell size is reduced. It would be identically zero 
once each cell is pierced by only one stellar line of sight, although the number density errors 
would grow as structure is washed out along the line of sight. 

3.3. Sampling density 

Finally we examined the ability of the algorithm to detect ISM structures for a limited 
stellar sampling density. Again we use the simulated ISM of Figure 5. We ran twelve mer 
reconstructions for N* from 13 to 342, with optimized cell numbers from 22 to 342 (Figure 7; 
Table 3). Using a column tolerance of 10~ 4 , the cell size was varied in increments of 5 pc to 
find the best reconstruction. Again, dark gray cells toward the edges which contain no stars 
are disregarded in the mer. 

The reconstruction using 13 stars can detect only the largest of the three clouds. The 
map is, however, quite good, considering that it gets the amplitude right, and matches the 
columns to within 0.003%. The reconstruction using 34 stars probably resolves the small 
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cloud closest to the sun. The 48-star map may actually resolve both small clouds, although 
it is not a robust detection. Although the column fit is nearly perfect, the map shows an 
unfortunate feature common to low-density sampled maps - the MER algorithm awards the 
pixel containing the sun an excess of material, at the expense of the 8 pixels which surround 
it. Thus if the goal of the reconstruction is to study the very local ISM, one cannot trust a 
high density contrast in only the central (solar) bin. By the time we get to 120 and 177 stars, 
both small clouds are clearly detected, as is the low-density bridge between the two largest 
clouds. Clearly, the sample of stars must be sufficiently dense such that the ISM structures 
to be detected are at least as large as the mean spacing between stars. 



4. Summary and applications 

Maximum entropy reconstruction is a powerful technique for recovering the structure 
of the ISM from column density data sets. By modulating the cell size, one can find a 
solution which recovers the underlying structure in the ISM. The lower limit to the size scale 
of detected structures is set by the stellar sampling density. The reconstruction algorithm 
is insensitive to the particular form of the entropy employed, so it is most prudent to use 
computationally simple and rapidly converging forms such as the quadratic entropy (equation 
14). 

The next step is to apply this technique to existing astronomical data sets. Different 
components of the ISM can be traced by absorption and emission measurments made at 
different wavelengths. Using a variety of absorptive and emissive tracers, we should be able 
to map regions of varying ionization and molecular state and elemental composition, and 
compare their spatial distributions. 

We are currently working to map the distribution of dust in the ISM using column 
densities derived from stellar color excess measurements. In addition, we are using archived 
ROSAT observations of hot stars to determine X-ray absorption columns. Since the soft X- 
ray cross section is insensitive to the molecular and (modest) ionization state of the absorber 
(Arabadjis & Bregman 1999), reconstructing the density field using X-ray columns will 
map the total density of material the local ISM. Additionally, by exploiting the spectral 
resolution of instruments aboard XMM-Newton and Chandra, one may be able to separate 
this material by ionization state. For example, AE/E ~ 50 for the non-dispersive ACIS 
detector aboard Chandra, resulting in a spectral resolution of about 10 eV at the prominent 
absorption edge of oxygen at 0.5 keV. Since the neutral and first three ionization states of 
oxygen are separated by about 20 eV each (the potentials are 13.6, 35.1, and 54.9 eV), it 
may be possible to distill the oxygen distribution into its dominant ionization states. 
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Several ISM tracers also contain velocity information. Optical absorption lines in stellar 
spectra, originating from the ISM, provide a measure of the distribution of metals not bound 
up in dust grains, and maps constructed from their columns may provide insight into the 
gas phase chemistry of ism. Additionally, the velocity structure of a resolved absorption line 
can be used in conjuction with a rotation curve and a second application of the maximum 
entropy principle to further constrain the distribution of material along the line of sight. 
Similarly, velocity-resolved 21 cm emission maps can be used to construct a distribution of 
neutral hydrogen along each line of sight. 

The applicability of this technique to absorption measurements at almost any wavelength 
comes at a particularly fortuitous time. Most astronomical satellite programs support an 
on-line archive containing their public domain observations, many of which can be mined for 
projected density measurements. With the number number of of these archives growing on 
an almost daily basis, this technique should find wide application in the near future. 

The authors would like to acknowledge support from NASA grant NAG5-3247. JSA 
would like to thank Doug Richstone for useful discussions. 
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Fig. 1. — Two discrete-bin ISM models which result in the same intervening column density 
toward a star. 




Fig. 2. — The four entropy forms examined in this study. Note the concavity inflection of 
the pseudo-entropy S = ne~ n at n — 2. 
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Fig. 3. — Simulated density map used to gauge the effect of changing the form of the entropy 
used in the reconstruction. The filled black dots represent the location of the stars whose 
intervening columns are used in each map reconstruction. The observer location is shown 
by the sun symbol at the map center. 
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Fig. 4. — Maximum entropy reconstructions of the density field shown in Figure 3 using four 
different entropy functions. The maps are nearly indistinguishable. 
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Fig. 5. — Simulated density field used to gauge the effect of modulating the cell size for a 
given stellar sampling density (Figure 6) and for modulating the stellar sampling density 
using an appropriate cell size (Figure 7). 




Fig. 6. — Maximum entropy reconstructions of the density field shown in Figure 5 using a 
range of cell sizes. 
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N* = 69 Nc = 124 N* = 89 Nc = 196 N* = 120 Nc = 233 




N* = 177 Nc = 358 N* = 272 Nc = 484 N* = 342 Nc = 489 




Fig. 7. — Maximum entropy reconstructions of the density field shown in Figure 5 using a 
range of stellar sampling densities and roughly optimized for cell size. 
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Table 1. Maximum entropy reconstructions of the simulated data set shown in Figure 3, for 
four different entropy functions. Relative errors are shown in parentheses. 



entropy 
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Table 2. Maximum (quadratic) entropy reconstructions of the simulated data set shown 
in Figure 5, for variable cell size s. Relative errors are shown in parentheses. 
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Table 3. Maximum (quadratic) entropy reconstructions of the simulated data set shown in 
Figure 5, for variable and optimized cell size. Relative errors are shown in parentheses. 
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0.0000 x 10 17 
(0.000 000) 


0.0000 x 10 18 
(0.000 000) 


(■ 


-2.540 
-0.1003) 


12.80 
(0.4770) 


34 


120 


43 


0.0000 x 10 17 
(0.000 000) 


0.0000 x 10 18 
(0.000 000) 


(- 


-0.5711 
-0.01688) 


8.198 
(0.2676) 


48 


80 


95 


0.0000 x 10 17 
(0.000 000) 


0.0000 x 10 18 
(0.000 000) 


(■ 


-2.826 
-0.1266) 


10.88 
(0.3787) 


69 


70 


124 


0.0000 x 10 17 
(0.000 000) 


0.0000 x 10 18 
(0.000 000) 


(■ 


-2.687 
-0.1169) 


10.74 
(0.3426) 


89 


55 


196 


0.0000 x 10 17 
(0.000 000) 


0.0000 x 10 18 
(0.000 000) 


( 


-2.544 
-1381) 


9.424 
(0.3753) 


120 


50 


233 


0.05 x 10 17 
(0.000 001) 


0.0258 x 10 18 
(0.000 003) 


(- 


-1.980 
-0.09216) 


8.327 
(0.3169) 


177 


40 


358 


0.011 x 10 17 
(0.000 000) 


0.0106 x 10 18 
(0.000 001) 


(- 


-1.692 
-0.08350) 


8.556 
(0.2961) 
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Table 3 — Continued 

N+ s N c (N mer -N) a/ ((N mer — N) 2 ) (n mer -n) V '((n mer - nf) 

(pc) (cm' 2 ) (cm' 2 ) (cm' 3 ) (cm' 3 ) 

272 35 484 2.217 x 10 17 1.907 x 10 18 -1.274 7.527 

(0.000 007) (0.000 029) (-0.05995) (0.2609) 

342 35 489 -1.447 x 10 17 18.60 x 10 18 -1.178 6.801 

(-0.000 041) (0.000 385) (-0.0597) (0.2549) 



